/*
- 
- Yuri Toledo Neves 8124361
- Trabalho de Recuperação SSC0143 - Programação Concorrente 
- Versão Paralela MPI
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "mpi.h"
#include <time.h>

//Funcoes de ajuste
double forcaArrasto (double H, double v);
double ajusteMassa (double t);
double ajusteForca (double t);
double ajusteGravidade (double H);

//EDOs 
double funcaoAceleracao(double t, double v, double gama, double H);
double funcaoRotacao(double t, double v, double gama, double H);
double funcaoHorizontal(double v, double gama);
double funcaoAltitude(double v, double gama);

//Método de Runge-Kutta de Ordem 4
void rungeKutta(int processId, int nprocess);

const double h = 0.0001;	// h do método de Runge-Kutta
const double g0 = 9.80665;	// Gravidade
const double R = 6335437.00000;	// Raio da terra
const double ro0 = 1.2922;	// Densidade do ar no nível do solo
const double hyp = -0.000101741; 	// Resultado equacao hypsométrica
const double PI = 3.141592;	// PI 
#define ROOT 0 			// Processo root

double	m0, mp, tq, fempuxo, alpha, coef_arrasto, sf, hf;	//Váriaveis de entrada
double 	r = 0.01; // Raio do foguete

//double gama;	// Angulo de lançamento	
double H = 0; 	// Altura 
double v = 0.1; // Velocidade
double t = 0;	// Tempo
double S = 0;	// Deslocamento vertical 
char op;
double erro;	// Erro

int main(int argc, char*argv[])
{	
	int i;
 	int rc,processId,nprocess;	//Váriaveis para controle de processos do MPI

	scanf("%c",&op);
	
		rc = MPI_Init(&argc, &argv);	// Verifica se os parametros necessários para inicialização do processo foram passados
		MPI_Comm_size(MPI_COMM_WORLD,&nprocess);	// Verifica a quantidade de processos
 		MPI_Comm_rank(MPI_COMM_WORLD, &processId);	// Enumera os processos
		MPI_Barrier(MPI_COMM_WORLD);				// Espera todos os processos chegaram a este ponto 
	
		if(processId == ROOT)
		{
			if(op == 'G')
			{	
				scanf("%lf %lf %lf %lf %lf %lf %lf %lf", &m0, &mp, &tq, &fempuxo, &alpha, &coef_arrasto, &r, &sf);
			}

			if(op == 'R')
			{	
				scanf("%lf %lf %lf %lf %lf %lf %lf", &m0, &mp, &tq, &fempuxo, &alpha, &coef_arrasto, &sf);
				r = 0.005;
			}

		}

		    MPI_Bcast(&m0, 1, MPI_DOUBLE, ROOT, MPI_COMM_WORLD);
		    MPI_Bcast(&mp, 1, MPI_DOUBLE, ROOT, MPI_COMM_WORLD);
		    MPI_Bcast(&tq, 1, MPI_DOUBLE, ROOT, MPI_COMM_WORLD);
		    MPI_Bcast(&fempuxo, 1, MPI_DOUBLE, ROOT, MPI_COMM_WORLD);
		    MPI_Bcast(&alpha, 1, MPI_DOUBLE, ROOT, MPI_COMM_WORLD);
		    MPI_Bcast(&coef_arrasto, 1, MPI_DOUBLE, ROOT, MPI_COMM_WORLD);
		    MPI_Bcast(&r, 1, MPI_DOUBLE, ROOT, MPI_COMM_WORLD);
		    MPI_Bcast(&sf, 1, MPI_DOUBLE, ROOT, MPI_COMM_WORLD);

		if(rc == MPI_SUCCESS)
		{
			
			MPI_Barrier(MPI_COMM_WORLD);  
			rungeKutta(processId, nprocess);
		}
			
		 MPI_Barrier(MPI_COMM_WORLD);
		 MPI_Finalize();
	return 0;
}

double forcaArrasto (double H, double v)
{

	double ro = ro0*exp(hyp*H);		// Densidade
	double D  = 0.5*(coef_arrasto*PI*r*r*ro*v*v);
	
	return D;

}

double ajusteMassa (double t)
{
	if(t < tq)
		return(m0 - mp*(t/tq));
	else
		return(m0 - mp);
}

double ajusteForca (double t)
{
	if(t <= tq)
		return(fempuxo);
	else
		return(0);
}

double ajusteGravidade (double H)
{
	return (R*R*g0)/((R+H)*(R+H));
}

//EDO da Aceleração
double funcaoAceleracao(double t, double v, double gama, double H)
{
	double arrasto = forcaArrasto(H, v);
	double massa = ajusteMassa(t);
	double aceleracao = ((ajusteForca(t)*cos(alpha))/massa) - (arrasto/massa) - (ajusteGravidade(H)*sin(gama));

	return aceleracao;
}

//EDO da Rotação
double funcaoRotacao(double t, double v, double gama, double H)
{
	double rotacao = (((ajusteForca(t)*sin(alpha))/(ajusteMassa(t)*v)) - ajusteGravidade(H))*(cos(gama)/v);
	return rotacao;
}

//EDO deslocamento Horizontal
double funcaoHorizontal(double v, double gama)
{
	double horizontal =  v*cos(gama);
	//printf("Gama %lf - Cosseno %lf\n",gama, cos(gama));
	return horizontal;
}

//EDO Altitude
double funcaoAltitude(double v, double gama)
{
	double altitude = v*sin(gama);
	//printf("Gama %lf - Seno %lf\n",gama, sin(gama));
	return altitude;
}

// Runge-Kutta de Ordem 4
void rungeKutta(int processId, int nprocess)
{
	double k1, k2, k3, k4; // Coeficientes do método para funcao Aceleracao
	double l1, l2, l3, l4; // Coeficientes do método para funcao Rotacao
	double m1, m2, m3, m4; // Coeficientes do método para funcao Horizontal
	double q1, q2, q3, q4; // Coeficientes do método para funcao Altitude 
	int i,j, flag = 0, found = 0;
	int num = processId;
	double g;	
	double g_graus;
	double altura = 0;
	double v1 = v, t1 = t, S1 = S, H1 = H;	// Váriaveis para proteção dos valores das variaveis globais
	float begin, end; 	// Determina o intervalo dos angulos
	double raio = r;
	begin =  num * (180/nprocess);
	end = (num+1)*(180/nprocess);
	g_graus = begin;

	do{
		
		if(op == 'R')
		{	
			if(g_graus > 70)
			{
				g_graus = 65;
				raio += 0.01;
			}
		}

		altura = 0;
		v1 = 0.1; // Velocidade
		t1 = 0;	// Tempo
		S1 = 0;	// Deslocamento vertical 
		H1 = 0;
		g = (g_graus*PI)/180;	// Converte o angulo de graus para radianos 
		
		for(i = 0; altura >= 0 ; i++)
		{
			
			k1 = funcaoAceleracao(t1, v1, g, H1);
			l1 = funcaoRotacao(t1, v1, g, H1);
			m1 = funcaoHorizontal(v1, g);
			q1 = funcaoAltitude(v1, g);

			k2 = funcaoAceleracao(t1+(h/2), v1+(h/2)*k1, g+(h/2)*l1, H1+(h/2)*q1);
			l2 = funcaoRotacao(t1+(h/2), v1+(h/2)*k1, g+(h/2)*l1, H1+(h/2)*q1);
			m2 = funcaoHorizontal(v1+(h/2)*k1, g+(h/2)*l1);
			q2 = funcaoAltitude(v1+(h/2)*k1, g+(h/2)*l1);

			k3 = funcaoAceleracao(t1+(h/2), v1+(h/2)*k2, g+(h/2)*l2, H1+(h/2)*q2);
			l3 = funcaoRotacao(t1+(h/2), v1+(h/2)*k2, g+(h/2)*l2, H1+(h/2)*q2);
			m3 = funcaoHorizontal(v1+(h/2)*k2, g+(h/2)*l2);
			q3 = funcaoAltitude(v1+(h/2)*k2, g+(h/2)*l2);

			k4 = funcaoAceleracao(t1+h, v1+(h*k3), g+(h*l3), H1+(h*q3));
			l4 = funcaoRotacao(t1+h, v1+(h*k3), g+(h*l3), H1+(h*q3));
			m4 = funcaoHorizontal(v1+(h*k3), g+(h*l3));
			q4 = funcaoAltitude(v1+(h*k3), g+(h*l3));

			t1 = t1 + h;
			v1 = v1 + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
			g = g + (h/6)*(l1 + 2*l2 + 2*l3 + l4);
			S1 = S1 + (h/6)*(m1 + 2*m2 + 2*m3 + m4);
			H1 = H1 + (h/6)*(q1 + 2*q2 + 2*q3 + q4);

			altura = H1;

		}
		
		g_graus += 1;

		if(g_graus >= end)
		{	
			flag = 1;
		}

	}while((fabs((S1 - sf)/sf) > 0.01) && (g_graus < end));

	if(flag == 0)
	{
		erro = fabs((S1 - sf)/sf);
		printf("%c %.5lf %.5lf\n",op, (g_graus-1)*PI/180, r);
		MPI_Abort(MPI_COMM_WORLD, 99);
	}

}
	